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elsewhere when an explicit method is used, are avoided. Moreover the method (an event- capturing, 
or time-stepping algorithm) allows for accumulation of events (Zeno phenomena) and for multiple 
switching surfaces (i.e., a sHding surface of codimension ^ 2). The details of the implementation 
are given, and numerical examples illustrate the developments. This method may be an alternative 
method for chattering suppression, keeping the intrinsic discontinuous nature of the dynamics on 
the sliding surfaces. Links with discrete-time sliding mode controllers are studied. 
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Simulations numeriques par la methode d'Euler implicite des 

systemes a modes glissants 

Resume : Dans ce rapport, on montre que la discretisation en temps de type Euler impli- 
cite conduit a une tres bonne stabilisation d'une classe de systemes commutes avec des modes 
glissants, et de leurs derivees sur la surface de glissement. Les oscillations artificielles qui sont 
generalement mentionnees pour I'implementation discrete de ce type de systemes sont evitees. De 
plus, la methode (de type "event-capturing" ou "time-stepping") permet de traiter des accumula- 
tions d'evenements (Phenomene de Zenon) et des surfaces de commutations multiples {i.e. des 
surfaces de glissement de codimension > 2). Dans ce rapport, les details de I'implementation sont 
donnes et des exemples numeriques illustrent ses proprietes. Cette methode pent (3tre une alterna- 
tive aux methodes complexes de suppression des oscillations, en gardant la nature intrinsequement 
discontinue de la dynamique sur les surfaces de gHssement. Le lien avec les commandes a modes 
glissants en temps discret est etudie. 

Mots-cles : Systemes commutes. Inclusion Differentielles de Filippov, Problemes de complemen- 
tarite, Methode d'Euler implicite, modes glissants, operateurs, maximaux monotones, Probleme 
lineaire de complementarite mixte. Discretisation Bloqueur d'Ordre Zero (BOZ) 
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1 Introduction 

Sliding mode controllers are widely used because of their intrinsic robustness properties [4T| [23t 
m [51] . Some important fields of application are induction motors |43[ [53t [7], aircraft control 
[44 } [32 j [54 } 135). hard disk drives |33[ I31j. solar energy systems [28]. However they are known to 
generate chattering which renders their appHcation deHcate. Solutions to cope with chattering 
or reduce its effects have been proposed, see e.g. [U O [l2l [H EH [55], which also have their 
own limitations [55]. One drawback of these solutions is that they usually destroy the intrinsic 
discontinuous nature of sliding mode control. Fundamentally, these control schemes are of the 
switching discontinuous type and they yield closed-loop systems that can be recast into Filippov's 
differential inclusions. The numerical simulation of such nonsmooth dynamical systems is non 
trivial and it has received a lot of attention, see e.g. [50l [HI [161 [HH [HI [371 [iQ], to cite a few. 
Both event-driven methods and time-stepping methods have been developed, see e.g. [l] for a 
survey. In this paper we focus on time-stepping methods, which have an interest not only for the 
sake of numerical simulation, but also for the real implementations of sliding mode controllers on 
discrete-time systems [52]. Recently it has been shown that the explicit Euler method generates 
unwanted effects like spurious oscillations (also called chattering effects) around the switching 
surface [251 [Ml [IS [57] . In parallel, the digital implementation of sliding mode controllers has been 
studied thoroughly in [27l I36j. where the Zero-Order Holder (ZOH) discretization is used. 

The purpose of this paper is to analyze the implicit {backward) Euler method for some particular 
classes of differential inclusions, that include sliding mode controllers. It is shown that, besides 
convergence and order results, the advantage of the impHcit method is that it allows one to get a 
very accurate and smooth stabilization on the switching surface (of codimension one or larger than 
one). Roughly speaking, this is due to the fact that the switches are no longer monitored by the 
state at step k, but by a multiplier (a slack variable in a nonlinear programming language). The 
multivalued part of the sgn(-) function, i.e. a multifunction, is then correctly taken into account, 
avoiding stiff problems. The advantage of such "dual" methods in terms of their accuracy on the 
sHding surface has already been noticed in [l9l[50] in an event-driven context, where the motivation 
was the simulation of mechanical systems with Coulomb friction. From a numerical point of view, 
our study shows that convergence and order results may not be sufficient to guarantee that the 
derivative of the state is correctly approximated on the switching surface. The implicit method 
adapts naturally to an arbitrary large number of switching surfaces, that is not the case of most of 
the other methods which become quite cumbersome as soon as more than two switching surfaces 
are considered. A further advantage of the proposed method is that contrary to other methods that 
have been studied and which destroy the intrinsic discontinuous nature of sliding mode systems 
(Hke the so-called boundary layer control, or various filtering techniques), our method keeps 
the multivalued discontinuity and consequently the fundamental aspects and properties of sliding 
mode control from a Filippov's systems point of view. Moreover, sampling rates need not be high 
to reduce chattering, contrary to other discrete sliding mode controllers. A second contribution 
of this paper is to show that the results that hold for the backward Euler scheme, extend to ZOH 
discretizations of sliding mode systems. 

The paper is organized as follows: Section[2]presents a motivating example for using an implicit 
Euler implementation of the simplest sHding mode system. In Section[3l a class of differential inclu- 
sions is introduced and existence and uniqueness results are given under the maximal monotonic- 
ity assumption. Through several examples, the Equivalent-Control-Based Sliding-Mode-Control 
(ECB-SMC) and the Lyapunov-based discontinuous robust control are shown to fit well within 
this class of differential inclusion. In Section [H some convergence and chattering free finite-time 
stabilization results are given. These central results of the paper show that the implicit Euler 
implementation of the differential inclusion yields a chattering free convergence in finite time on 
the sliding surface. Section [5] is devoted the study of Discrete-time SHding Mode Control and the 
extension to ZOH discretization. Some hints on the numerical implementation of the implicit Euler 
scheme are given in Section [6] and the paper ends with some numerical experiments in Section [9l 

^see [55] for a discussion on this point. 
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Notations and definitions: Let A G IR^^™', then A,i is the ith column and Ai, is the ith 
row. The open ball of radius r > centered at a point x € iR" is denoted by Brix). For a set 
of indices a C {1, . . . ,rt} and a column vector x G IR^, the column vector Xa will denoted the 
sub- vector of corresponding indices in a, that is Xa = [xi,i G a]'^ . 



2 A simple example 

To start with we consider the simplest case: 

r 1 if x{t) < 
x{t) e -sgn(x(i)) = < -1 if xlt) > , a;(0) = xq (1) 
[ [-1,1] if x{t)=0 

with x{t) £ IR. This system possesses a unique Lipschitz continuous solution for any xq- The 
backward Euler discretization of JT]) reads as: 

Xk+i - Xk ^ -hsk+i 
jk+i e sgn(xfe+i) 

This method converges with at least order ^ (see Proposition [2] below) . Let us now state a result 
which shows that once the iterate Xk has reached a value inside some threshold around zero for 
some k, then the dual variable Sk+i keeps its value and so does Xk+n for all n ^ 1. 

Lemma 1 For all h > and xq G M, there exists ko such that Xko+n — and ^"+"+^ fco+n _ 



for all n ^ 1. 

Proof: The value fco is defined as the first time step such that Xko & [—h, h]. If xq e [—ft,, h], 
then kf) = 0. Otherwise, the solution of the time-discretization ^ is given by Xk = xo — 
sgn{xo)kh, Sk = sgn(xo) while Xk ^ [—h,h] for k < ko, and fco = Tt^I ~ 1- The symbol 
[x] is the ceiling function which gives the smallest integer greater than or equal to x. 

Let us now consider that Xkg G [—h, h]. The only possible solution for 

Xko + l ~ Xko = ^hsk„ + i 

(3) 

Sfco+i e sgn{xk,,+i) 



is Xko+i ~ and Skg+i ~ For the next iteration, we have to solve 



Xko+2 = ~hSko+2 

Sko+2 e sgn(a;feo+2) 

and we obtain Xko+2 ~ and Ska+2 = 0. The same holds for all Xka+mSka+n, ri ^ 3, redoing the 

same reasoning. Clearly then the terms ^^0+"+^ — Xk^+n ^ppj-Qxiniating the derivative, are zero 

h 

for any h > 0. □ 
This result is robust with respect to the numerical threshold that can be encountered in floating 
point operations. Indeed, let us assume that Xkg — h = e ^ 1, that is, e > is zero at the machine 
precision. We obtain Sk„+i = — 1 and Xk^+i = £ that is zero at the machine precision. For n = 2, 
we obtain Xko+2 = and Sko = J^- This robustness stems from the fact that the dynamics is not 
only monitored by the sign of Xk but also by the fact that the "dual" variable Sk+i belongs to 
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Figure 1: Iterations of the backward Euler method. 



Consequently this result shows that there are no spurious oscillations around the switching 
surface, contrary to other time-stepping schemes Hke the explicit Euler method [HI [26]. Remark- 
ably Lemma [U holds for any > 0, which means that even a large time step assures a smooth 
stabilization on the sliding surface. It is noteworthy that solving the system ([2]) with unknown 
Xk+i and Sfc+i is equivalent to calculate the intersection between the graph of the multivalued 
mapping Xk+i ^ — /isgn(xfe+i) and the straight line Xk+i ^ Xk+i — Xk- This is illustrated on 
Figure [U where few iterations are depicted until the state reaches zero. 

Prom a control perspective the input is implemented on [tk,tk+i) as = — sgn(xfe+i) as a 
function of Xk and /i, where h is the sampling time. There is no problem of causality in such 
an implementation. It is noteworthy that in the implicit method there is absolutely no issue 
related to calculating sgn(O), or more exactly sgn(e) where e is a very small quantity whose sign is 
uncertain. The implicit method automatically computes a value inside the multivalued part of the 
sign multifunction and may be considered as the time-discretization of the multifunction sgn(-). 
It is easy to show that the explicit method yields an oscillation around a; = 0, as shown in more 
general situations in |251I26). Other time-stepping methods exist, like the so-called switched model 
[U [37] , however it fails to correctly solve the integration problem when the number of switched 
surfaces is too large (see also [1] for similar issues when the so-called sigmoid blending mechanism 
is implemented). Moreover this method may yield a stiff system, and from a control point of view 
it introduces a high-gain feedback that may not be desirable in practical applications. 

On Figure[2]Ia)-(c), the discrete state x^ and the control are displayed for xq = 1.01 at to = 
and for various values of the time-step h that are sufficiently large to illustrate the behavior of 
the time-stepping scheme and its convergence. 

Let us define two discrete function norms to measure the convergence: 

l|e/lloo = Ej=o 1//^ - fi*k)\ 

(5) 

\\ef\\p = ihj:lo\h~fit'^)n^''- 

We can compute that 

||e,|U = 1 for all > (6) 
and therefore there is no convergence in infinite norm ||.||oo for s = sgn(a;). In ||.||i and ||.||2, we 



can respectively observe the convergence with order 1 on Figure 2(d) 



Complementarity framework Let us end this section by restating the systems H]) and ^ 
into the complementarity framework. Let us introduce equivalent formulations of the inclusion 
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Figure 2: A simple example for xq = 1.01 at to = 0. 
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s{t) e sgn(a;(i)) such that 

( x{t) = \{s{t) e] - l,i[ 

s(t) e sgn{x{t)) ^ x{t) e iV[„i,i](s(t)) s{t) e [-1, l] and I x{t) if x{t) = -1 (7) 

\x{t) ^ iix(t) = 1 

where -/V[_i.i] is the normal cone in the sense of Convex Analysis to the interval [—1,1]. The 
definition of the normal cone in the present case, 

= {-vi + V2, i>i _L s + 1 ^ 0, W2 -L 1 - s ^ 0} (8) 

yields the following complementarity representation of the sign multi- valued function 

{x{t) ^ ~Vi{t) ^ V2{t) 

x{t) e iV[_i,i](s(t)) 4^ < s$ -yi (t) _L s{t) + 1^0 (9) 
[O S$ V2{t) -L 1 - s{t) ^ 

In order to directly substitute the value of s{t) into the dynamics x{t) = —s{t), a other comple- 
mentarity formulation can be defined. By setting Xi{t) = ^(1 — s(t)) and X2{t) = vi{t), one gets 

rs(t) = l-2Ai(t) 

x{t) e ^ <^ ^ Ai(t) _L x{t) + A2(t) ^ (10) 

[o ^ X2{t) ± 1 - Xi{t) 
3 A class of differential inclusions 

Let us now introduce the following class of differential inclusions, where x{t) G M": 

i{t) g -A{x{t)) + f{t,x{t)), a.e. on (0,T) 

(11) 

a;(0) = xo 
The following assumption is made: 

Assumption 1 The following hold: 

• (i) A{-) is a multivalued maximal monotone operator from J?" into IR^ , with domain D{A), 
i.e., for all x £ D{A),y G D{A) and all x' G A{x),y' G A{y), one has 

{x' -v')^{x-y)^Q (12) 

• (ii) There exists L ^ such that for all t G [0,r], for all xi,X2 G 1R'\ one has \\f{t,xi) — 

fit,X2)\\^L\\xi-X2\\. 

• (iii) There exists a function <!>(•) such that for all R ^ 0; 

$(i?) = sup jll ^i-^v) \\cmo,Ty,lR") I II ^ hmo, T);1R")^ ^} < +°° 

The following is proved in [TOl [9]. 

Proposition 1 Let AssumptionU\hold, and let xq G D{A). Then the differential inclusion (Qlp 
has a unique solution x : (0, T) J?" that is Lipschitz continuous. 
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In this paper we shall focus on inclusions of the form: 

g f{t, x{t)) - BSgn{Cx{t) + D), a.e. on (0, T) 

x(0) = 

A 



(13) 



with B e K"''", and Sgn{Cx{t) + D) = (sgn(Cia; + Di), sgn(C,„a; + Dm)f G -K". It will be 
shown how to recast lfT3|) into (fTTj) . 

Example 1 (Equivalent-control-based sliding-mode-control (ECB-SMC)) Consider a sys- 
tem i(t) — Fx{t) + Gu, with an equivalent- control-based sliding-mode-control (ECB-SMC) of the 
form u{x) = -{HG)-^HFx - a{HG)-^ Sgn{Hx) , a > (see e.g. J57^). Then the closed-loop 
system x{t) = (F - G{HG)-^HF)x(t) - aG{HG)-'^ Sgn{Hx(t)) fits within 

Let us now state a well-posedness result which is a consequence of Proposition [H 

Corollary 1 Consider the differential inclusion in fl3]) . Suppose that (ii) and (iii) of Assumption 
QP hold. If there exists an n x n matrix P = > such that 



PB,^ = Gl 



(14) 



for all 1 ^ i ^ m, then for any initial data the differential inclusion / Qgj) has a unique solution 
X : (0, T) that is Lipschitz continuous. 

Proof: The proof uses a state variable change introduced in [13]. Let R be the symmetric 
square root of P, i.e. = P. Let us perform the state transformation z = Rx. Then we get 



ilt) e Rf{t, R-^z{t)) - RBSgn{GR-^z{t) + D) 



(15) 



Notice that BSgn{G R-'^ z{t) + D) ^ J2Zi B,.,sgn{G,R-^ z + A). Therefore RBSgn{GR-^z{t) + 
D) = E™i-R-B.,sgn(Q.i?-iz + A) = ET=iB'^ClsgniG,^R-^z + A)- We can rewrite the 

system as 



z{t) e Rf{t, R-^z{t)) - J2 R~'GZsgniG,,R-^zit) + A) 



(16) 



The multivalued mapping ^ i-^ sgn(^) is monotone. By [l6l Exercise 12.4] it follows that each 
multivalued mapping z i-^ R^^Gj,sgii{Gi,R~^z{t) + Di) is monotone. From [291 Proposition 
1.3.11] it follows that R~^Glsgn{G„R~^ z{t) + A) = dfi{z) with f,{z) = \Gi,R-^z{t) + Di\. By 
[45t Theorem 5.7] it follows that fi{-) is convex. Being the subdifferential of a convex function, 
the multivalued mapping z i-^ dfi{z) is maximal (monotone) [HI Corollary 31.5.2]. Therefore by 
Proposition [T] the inclusion in (fT6l) possesses a unique Lipschitz solution on (0,T) for any T > 
and since R is full-rank so does ifTSl). □ 



Example 2 Consider the 



mode system in f25\ Equ.(l)-(4)]. One has B = {0 1)-^, 



C = (ci 1), D ^ 0. Then the condition in holds with P = 
assures that P > 0. 

2 



Pll Cl 

ci 1 



and pii > (ci)^ 



Example 3 Consider B 



, Sgn{Gx + D) 



1 

2 -1 

tories may slide on codimension one surfaces xi + 2x2 ~ or 2xi 
2 surface (xi + 2x2 = and 2xi — X2 = 0). 



{sgn{xi + 2x2), sgn{2xi — X2)Y' . Trajec- 
X2 ~ and on the codimension 



Example 4 One solution to reduce chattering is the observer based SMC. Let us consider the 
following example taken from J55^ . whose closed-loop dynamics is given by: 



/ x{t) \ ( 
e{t) 

Xs{t) 














k 


-k 
















1 


V ^ 





_ 1 


_ 2 
r 



\ 


/ 


x{t) \ 










e{t) 











Xs{t) 







/ 




Xs{t) ) 




I / 



sgn{Gx{t)) 



(17) 
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with C = (1 — 1 0). For the notations see f55i §11.0]. This system satisfies the condition 
/ 1 -1 \ 



with P 



-1 P22 
P33 

V Pii J 



, P22 > 1, P33 > 0, P22 > 0. 



Notice that the condition implies that Bl^PB,i = BlfiJ, = Bi,C,i > 0. When m ^ 1 
this is a relative degree one condition. It is noteworthy that l|14p does not imply that B has full 
column rank. In particular it does not preclude m > n. Dissipative systems with no feedthrough 
matrix satisfy an input-output constraint similar to (fT4|) [T4] . 

Example 5 (Lyapunov-based discontinuous robust control) Let us show how the above ma- 
terial adapts to this type of feedback controller. The class of dynamical systems is 

xit) = f{x{t)) + Bu{t) + B-f{t), x{Q) = a;o (18) 

where x{t) € iR", B € IR^^™ , /(•) satisfies assumptionUl andj{-) S 5?™ is a bounded disturbance 
satisfying \ji{t)\ < pi for all I ^ i ^ m, all t ^ and some finite pi. The problem is the 
stabilization of the system at the origin x — 0, knowing that there exists a function V{-) such 
that the uncontrolled undisturbed system x{t) ~ f(x{t)) admits V{-) as a Lyapunov function. In 
particular, one has V{x{t)) — \7V{x{t))'^ f{x{t)) ^ along the trajectories of the free system. Let 
us rewrite the system in ^8\) as 

m m 

X{t) = f{x{t)) + B,^U, + B,^J^{t) (19) 

1=1 1=1 

Let us propose the control input Ui{x) = — pisgniVV'^ {x)B,i) . We obtain: 

ni m 

x{t) e f{x{t)) - Y P^B,,sgn{WV{xfB,,) + ^ S.,7,(t) (20) 

i=l i=l 

We can state the following result. 

Corollary 2 Suppose that V{x) = ^x^ Px, P = P^ > 0. The system in f20\) has a unique 
Lipschitz continuous solution on [0, +oo) for any xq. 

Proof: We have WV{x)'^ B,i = Bj^Px. Let z = Rx, where R > is the symmetric square 
root of P. We may rewrite f20]) as 



m e Rf{R'^z{t)) -J2p^RB,,sgn{BlRz) + ^i?B.,7,(t) 

j=i j=i 

Then following the same steps as for the proof of Corollary]^ we conclude that Proposition]^ applies 
to this system, hence to 1120]} . □ 
Such a controller assures the global asymptotic stability of the equilibrium x = 0. This is 
made possible because of the multivalued characteristic of the discontinuous input. The closed-loop 
system possesses the origin as its unique equilibrium, because of the multivaluedness property. The 
restriction to quadratic Lyapunov functions stems from monotonicity preserving conditions, and 
is not straightforwardly avoided. 

4 Convergence results and Chattering Free Finite-time Sta- 
bilization 



The differential inclusion ifTTj) is time-discretized on [0, T] with a backward Euler scheme as follows: 

(21) 



+ Aixk+i) B fitk,x,), for all k € {0, ...,N - 1} 
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where h= ^. The fully implicit method uses f{tk+i,Xk+i) instead of f{tk,Xk)- The convergence 
and order results stated in Proposition [2] below have been derived for the semi-implicit scheme 
(f2T1) in [l^. So the analysis in this section is based on such a discretization. However this is only 
a particular case of a more general 0— method which is used in practical implementations. The 
next result is proved in [TO] . 

Proposition 2 Under assumption[^ there exists rj such that for all h > one has 

For allt e[0,T], \\x(t) ~ it)\\ Vh (22) 
Moreover limft^o+ maxtg[o,T] \\xit) - + \\x{s) - x^{s)\\'^ds = 0. 

Thus the numerical scheme in l(2T|) has at least order and convergence holds. The conditions 
of Assumption [T] for the existence and uniqueness results of Proposition [T] are sufficient only. Other 
criteria exist, like Filippov's criterion for uniqueness of solutions [iHl Proposition 5]. Similarly it 
is possible that time-stepping methods converge for systems that satisfy such a criterion, despite 
no result seems to be available in the literature. As seen in Lemma[Tl the precision of the method 
may be much better than what is to be expected from ((22l) on large portions of the trajectories. 

The differential inclusion in (flSl) is therefore discretized as follows: 



( Xk+i Xk ^ _^(^^^^^) _ BSgii{Cxk+i + D), a.e. on (0,r) 

(23) 

^ x(0) = Xq 

One sees that advancing the implicit method from step k to step fc -I- 1 involves solving generalized 
equations with unknown Xk+i, of the form £ Fs{xk+i) + Fm{xk+i) where Fs{-) is singlevalued 
while Fm{-) is multivalued, h, tk and Xk appear as parameters of the generalized equations. Solving 
such generalized equations thus boils down to computing the intersection between the graph of 
Fs{-) and the graph of fm(-) as illustrated in section [H The result of Proposition [2] applies to 
((23|) . As we shall see next, such an implicit method also assures a good estimate of the derivative 
X and a smooth stabilization of the discrete-time solution on the sliding surface. 

Before stating the smooth stabilization result, let us consider a preliminary result. Let us 
denote the output of the dynamical as: 

y{t) ^ Cx{t) + D (24) 



Lemma 2 Let us assume that a sliding mode exists for some indices i S a C {1 . . .to} such that 

3<, > 0, ya{t) = Ca,x{t)^Da = 0, forallt>t^. (25) 

Then there exists p > such that for all t > t,, and for all x{t) such that Ca,x{t) + Da = 0, one 
has 

\\{Ca,f{x{t),t))\\^p (26) 

Furthermore, let Assumption l.(ii) holds, then the following hound is satisfied in the neighborhood 
of the sliding subspace, 

3r > 0,3k > 0,3p > such thafit > t*,Vx G Br{x), \\{Ca,f{x,t))\\ Kr + p (27) 

for all x(t) such that Ca,x{t) + = 0. 

Proof: From we have ya{t) e Ca,f{x{t),t) - hCa,BSgn{y{t)). For t > t*, the sliding 

mode ya{t) = implies that ya{t) = Ca,x{t) = for all t > t* and therefore 

Ca.fix{t), t) e C„.SSgn(y(t)) (28) 
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The inclusion |[28]) yields 

3p>0,|l(C„./(x(t),t))|Kp (29) 

for all x{t) such that Ca,x{t) + Da = 0. By the assumption l.(ii), the Lipschitz continuity of 
/(•,•) allows us to write for some k > 

V.f (i) e Br{x{t)), \\Ca,{f{^{t),t) - f{x{t),t))\\ < \\Ca,\\Lr = KT. (30) 

Combining ((29)) and ((3Q|) ends the proof. □ 
Lemma [T] extends to ((23)) as follows when the sliding surface of codimension \a\ is attained. 

Lemma 3 Let us assume that a sliding mode occurs for the index a C {1 . . - m}, that is ya{t) ~ 
0,t > t^. Let C and B he such that \14^ holds and Ca,B,a > 0. Then there exists he > such 
that V/i < he, there exists fco G iV such that Uko+n = Cxk„+n+i + D = for all integers n ^ I. 

Proof: At each time-step, we have to solve for yk+i — Cxk+i + D and Sk+i the generaHzed 
equation 

f Uk+i = Vk + hCf{tk,Xk) - hCBsk+i 

\ (31) 
[ Sk+i e Sgn(?/fc+i) 

Under condition l(H| , the convergence of the time-stepping scheme is ensured by Proposition ID 
The convergence and the existence of the sliding mode ensure that 

3fco,3A'i > 0,37^2 > 0,3ti > such that \\yaM-,\\ «^ KiVh and \\xko ~ x{ti)\\ A'zV/I (32) 

for Ca,x{ti) + Da = 0. Using ((27)) for x{ti) and a sufficiently small h such that r = K2Vh, we 
have the following bound 

||ya,fco + hCa,,f{tko,Xko)\\ Vh{Ki + hKK2 + Vhp) (33) 

Introducing the complementary index set /3 = {i,yi{t) = Ci,x(t) + Di ^ 0}, for t > t^, almost 
everywhere and using ((33)) we obtain that there exists pi > such that 



WVa.ko + hCa,,fitko,Xko) - hCa,B,i3Sgn{yi3^ko+i)\\ =^ Vh{Ki + hKK2 + Vh{p + pi)) (34) 
and therefore it is possible to choose hi such that for all h < hi 

I [-h{Ca,B,a)~^ [ya,ko + hC a., f {tko, ^ko) ^ /iC'q, S,/3Sgn(y/3,fco+i )]] J 1, for all i G a. (35) 
If 1(35)) is satisfied, the unique solution of l(3T)) at the iteration /cq + 1 is given by 

ya.ko + l = 0; Sa,ko + l = -h{Ca,B,ay^ [Va,ko + hCa,, f [tka , Xko) - hCa,B,[3Sgn{yf3,ko + l)] (36) 

The next iterate will by given by the solution of the generalized equation, 

J yko+2 = hCf(tko+i,XkQ+i) - hCBsko+2 
\ Sfeo+2 e Sgn(?/fc„+2) 



(37) 



Using the fact that ya,ko+i = Ca»Xko+i + Da — 0, we can use ((29| to conclude that there exists 
/i2 such that for all h < h2 

I [-h{Ca,B,ay^ [hCa,,f{tko+iiXko+i) " hC a.B,pSgn[y p ^ko+2)\\ ^\ ^ 1, for all i G a, (38) 
and therefore the solution of l|37p is 

yaM+2 = 0; SaM+2 = -h{Ca,B,ay^ [kCa,, f {tko + 1 ^ Xko + l) - hCa.B,i3Sgn{yi3^ko+2)] (39) 



The bound ((29)) is uniform and can be applied for the next steps. Choosing he as the minimum 
of the considered time steps hi, /12, ■ • ■, the proof is obtained for ya,ko+ni ri ^ I. □ 
The finite-time convergence of the time-discretization of similar nonsmooth dynamical systems 
(essentially mechanical systems with dry friction) is proved in [6]. Our results may therefore 
be considered as the continuation of studies on the finite-time convergence for algorithms of the 
proximal type. 
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5 Discrete-time Sliding Mode Control (SMC) 

This section is devoted to show how the above discretizations may be used in a digital control 
framework. 



5.1 Example of an implicit Euler controller (lEC) 

Let us come back to the inclusion in For this simple system, the ZOH and the Euler discretiza- 
tion yield the same-discrete system. Assume the integrator x{t) = u{t) is sampled with sampling 
period h > 0. On the time interval [tk,tk+i) one has x{t) = Xk + hfUk, where ht = t — t^. The 
controller u{x) = — sgn(a;) is known as the equivalent control-based SMC [57]- Let us implement 
a "backward" controller = — sgn(a;fe+i) at time tk, following the above lines. Suppose that 
Xk e [—h, h]. Then following the same calculations as in the proof of Lemma [U we obtain that 
Sfc+i = X- Therefore on [tk,tk+i): 

x{t) =Xk- j-xk (40) 
and it follows that x{tk+i) = Xk+i = 0. On the next sampling interval [tk+i,tk+2) one obtains 

Sfe+2 = 

x(0=Xfc+i-^Xfe+i= 0-^0 = (41) 

and so on on the next intervals, where the zero value is obviously some small value at the machine 
accuracy, if we suppose that Xk ^ [—h, h], the value of Sk+i is 1 or —1 according to the sign of Xk- 
To summarize the control is given explicitly in terms of Xk and h by 

Ufc = -proj[_i^i]( — ) (42) 

where proj,^ denotes the Euclidean projection operator onto the set C. 

As alluded to above, such an "implicit" input is causal and can be computed at tk with the 
values of the state at tk by l(42|) . It requires at each step to solve a rather simple multivalued 
problem which a Mixed Linear Complementarity Problem (MLCP, see Section [6] below) . It is not 
of the high gain type. 

Remark 1 The fact that the function sgn{-) generates only binary values (+1 or ~1) does not 
hamper the above method to work. Indeed the implicit Euler method allows us to compute values 
of the sign multifunction inside its multivalued part at Xk = 0. 



5.2 Extension to ZOH discretized systems 

The ZOH discretization of linear time invariant systems x{t) = Fx{t) + Gu{t) with an ECB-SMC 
controller, u{x) ~ —{CG)~^{CFx + aSgn(Ca;)),a > results in a discrete-time system of the 
form: 

Xk+i = ^Xk - Tsk for all t e [kh, {k + l)h) (43) 
where /i > is the sampling period, and 

$ = cxp(F/i)- / cxp{FT)dTG{CG)-^CF (44) 
Jo 

r = / eyi^{FT)G{CG)-^dT (45) 
Jo 

with G £ J?"^™, C G J?™^", when a explicit Euler implementation of the control is performed 
[52l[56]. For an implicit Euler implementation, let us set 

r Uk^-{CGr\CFxk + Sk+i) 

{ (46) 
I Sk+i = Sgn(Cx/c+i), 
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which corresponds to the implicit discrete time version of the ECB-SMC controller. We therefore 
get on each sampling period: 

Xk+i = $a;fc - Tsk+i for all t e [kh, {k + l)h) (47) 

At each time-step, one has to solve 

Xk+1 = - Tsk+i 

Vk+i ^ Cxk+i + D . (48) 
Sk+i € Sgn(yA;+i) 

Inserting the first line of (08]) into the second line we obtain the following one-step system 

Vk+i = C^Xk +D- CTsk+i 

(49) 

Sk+i e Sgn(yA;+i) 

Comparing with the time-discretized systems in ((23| and l|3ip one sees that the term hCB is 
replaced in case of a ZOH method by the term CT. Provided the problem has a unique solution 
one can compute the controller in l(i6l) with the knowledge of Xk, h, F, G and C. We will see in 
the next Section how the computation can be carried out in practice. 

6 Implementation of Discrete-Time Systems 

Let us consider in this section the following discrete-time system: 

Xk+1 = Rxk + P- Ssk+i 

Vk+i = Cxk+i + D (50) 
Sk+i e Sgn(2/fc+i) 

where /c ^ is an integer, Xk the discrete state, yk the discrete output and Sk the discrete input. 
The discrete system l(50|) is a common representative for the discretization given by l(23|) . l(2T|) or 
(081) and the matrices i?GiR"^",S'eiR"^™ and the vector pel?" are determined by the chosen 
time-discretization method and detailed in Section [Ql The matrices C and D are given by their 
definition in (fT3)) . 



6.1 Mixed Linear Complementarity Problem (MLCP) 

The time-discretized system (|50)) appears to be a Mixed Linear Complementarity Problem (MLCP) 
that we have to solve at each time-step. Let us define what is a MLCP in its general form with 
bounds constraints as it has been proposed in [2T] : 

Definition 1 (MLCP) Given a matrix M e 1R"^^"\ a vector q G 5?'" and lower and upper 
bounds l,u ^ M , find z e JR™, w, t; e such that 



Mz + q ~ w — V 

I !^ Z 

{z^l)'^w^O 
{u- z fv = 



(51) 



where IR = 1RU {+oo, — oo}. 
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Note that the problem ([51]) imphes that 

-(A/z + g) (52) 

where the notation Nc{x) is used for the normal cone in the Convex Analysis sense to a convex 
set C at the point x. The box [I, u] C 5?™ is defined by the Cartesian product of the intervals 
[li,Ui\,i G {1, . . . , to}. The normal cone to a convex set is a standard instance of a multi-valued 
mapping [l^ . The relation l(52|) is equivalent to the MCP l(5T|) if we assume that w is the positive 
part of Mz + g, that is w = [Mz + q)^ = max{0, Mz + q)) and v is the negative part of Mz + q, 
that isv = {Mz + q^ = niax{0, -{Mz + q)). 

In order to state the problem (|50|) as a MLCP, the variable Xk+i is condensed into the second 
line such that 

f yfc+i = CRxk + Cp- CSsk+i + D 



\ Sk+i e Sgn(yfc+i) 
and the following variable and parameters are defined as follows 



(53) 



z = sk+i; yk+i = w-v 

M ^CS, q = - [CRxk + Cp + D) (54) 

li — —l,Ui = 1,2 = 1 ... TO. 

Finally, the problem l(50|) can be recast into a MLCP by observing that 

Sfc+i e Sgn(yfe+i) 

t 

yk+i e A^[_i,i]™(sfe+i) 

t (55) 
Vj.k+i = if Sj,fc+i e] - 1, 1[ 

Sk+1 e [-1, 1]'" and ^ yj,k+i < if Sj^k+i = -1 , j e {1, ■ • -,'71} 

,y],k+i ^ if Sj,fc+i = 1 



The MLCP (|5T|) is a well-known problem in the mathematical programming theory arising 
for instance from the Karush-Kuhn-Tucker optimality conditions of a quadratic program or from 
the primal/dual optimality conditions of a Hnear program. The MCLP enjoys a large number 
of numerical algorithms and several rehable solvers have been implemented. Several famihes 
of solvers may be cited: a) extensions of Lemke and principal pivotal techniques for LCP to 
MLCP [m [m dH [n] b) extensions of projection/splitting techniques for MLCP [24l[19] and c) 
semi-smooth Newton methods [40]. In this paper, the computations are done with the help the 
SiCONOS/NuMERiCS Open source Library [3] and/or the PATH solver [21]. The results of existence 
and uniqueness of solutions of (jSTj) are related to the properties of M (P-properties or coherent 
orientations of the associated affine map (normal map) for particular cases of bounds constraints) . 
Without entering into further details, we refer to [30^ [24] for the main results. The assumptions 
on the matrix M drives the choice of particular solvers that can be in polynomial time rather than 
standard exponential time for brute force enumerative solvers. 



6.2 Some Time— Discretization Methods 

In this Section, the formulation of the discrete-time system l(50|) is related to the continuous time 
system l|13p through a given discretization method. 

Explicit Euler discretization of /(•,•) Let us start with the explicit Euler discretization 
method of the term f{t, x{t)) as it has been given in ((2T|) . At each time step, the matrices in l(50| 
and in the MLCP l(5T|) can be identified as 

R^I,p^hf{tk,Xk),S = hB, M = hCB,q^~ihCf{tk,Xk) + Cxk+D) (56) 
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Discrete-time Plant 



MLCP solver 



Figure 3: Control system schema with imphcit Euler implementation. 

Let the assumptions of Corollary [T] be satisfied with B full-column rank {CB ~ B^PB > 0). This 
result ensures the existence and uniqueness of a solution of the MCLP. Furthermore, standard 
pivotal techniques such as Lemke's method or projection/splitting such as Projected Successive 
Over-Relation (PSOR) compute the solution. 

Implicit Euler and 6- method In a more general way, we can choose to time-discretize the 
term f{t,x{t)) by a implicit Euler scheme or a 6'-method. The main motivation for doing in this 
way is the higher accuracy and stability that we can obtain for such a numerical integration scheme 
(see [2] for an example of instability with the ExpHcit Euler method) . Let us consider first that the 
mapping /(•, •) is afRne, that is f{t,x{t)) = Fx{t)+g. The matrices in ([50]) and in the MLCP l(5T1) 
can be identified as 

f R={I - h6F)-^{I + h{l - 0)F),p ={I - h6F)-^g, S = h{I - hOFy'^B, 

\ M = hC{I - heF)-^B,q= -{{I - h9F)-^{I + h{l - 9)F)xk + {I - hOF)-^ g + D) ^^'^'^ 

for 9 e [0, 1]. For 9 = 0, the explicit Euler case is retrieved. For 9 = 1, the implicit Euler scheme 
is used to discretize /. If the mapping /(•, •) is nonlinear, a newton linearization can be invoked. 
In this case, the solution at each time step is sought as a Hmit of solutions of successive MLCPs. 
We refer to [2] for a detailed presentation of these developments. 

Zero— Order Holder (ZOH) method The ZOH discretization presented in Section [5?2l can be 
also formaHzed into the form (|50)) and then ifST]) with 

i? = $,p = 0,5* = F, il/ = Cr,fy = -(C$Xfe) (58) 

In practice, numerous methods are available to compute the ZOH discretization, i.e., $ and F 
which amounts to compute the matrix exponential and its time integral [39]. In this work, the 
numerical computation is performed using an explicit Runge-Kutta method with high order of 
accuracy and a numerical tolerance near the machine precision threshold. On the Figure O the 
control scheme is depicted showing that the controller is causal and computed form Xk- 

7 Two other classes of differential inclusions 

In this section, we introduce other classes of differential inclusions which extend ifTS)) . The second 
class of differential inclusions is: 

x{t) e fit, x{t)) - YZi (^«a;(i) + Bi)sgTi{C,x{t) + A), a.e. on (0, T) 

(59) 

x(0) = xo 

with A, e K"""", B, e d e 1?^''", A is a scalar. 
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The third class that we shall analyze is: 

xit)efit,xit))-gixit))Sgnihixit)), a.e. on(0,T) 

(60) 

a;(0) = xo 

where g : 5?" _ZR"^™ and h : WC -ZR™ are smooth functions, Sgn(ft.(x)) = [sgn(/ii(a;), sgn(/im( 

Corollary 3 Consider the differential inclusion in Suppose that (ii) anrf (iii) of Assumption 

QP ftoZd. Suppose that the multivalued mappings x (A^a; + Bi)sgn{CiX + Di), 1 ^ i ^ m, are 
hypomonotone. Then for any initial data the differential inclusion (53) has a unique solution 
X : (0, T) WC^ that is Lipschitz continuous. 

The proof is straightforward and is omitted. 

Example 6 The mapping IR M, x t-^ [x + l)sgn{x) is hypomonotone. Indeed {x + l)sgn{x) = 
|a:| + sgn{x) and x ^ \x\ + kx is monotone for any k ^ 1. 

{-kix + 1 if x^O 
— k2X — 1 if X ^ ; is 
[-1,1] if x = 

hypomonotone with constant k for any k ^ max(|fci |, |fc2|). Let fci — — fc2 = k. Then F{x) = 
{kx + l)sgn{x). The linearized Stribeck friction model (with multivalued part at zero tangential 
velocity) f38^ is hypomonotone. 

Let us state other cases where l(59|) fits within Proposition [TJ 

Lemma 4 Let Bi — aCj for some a > Q, Di = Q and Ker{Ci) C Ker{Ai). Then for any 
initial data the differential inclusion / 1 59)) has a unique solution x : (0, T) IR^ that is Lipschitz 
continuous.. 

Proof: First notice that the set-valued mapping x ^ aCfsgii{CiX + Di) is maximal monotone 
[i6t Exercise 12.4]. Under the lemma's conditions, one sees that x i-^ AiXSgn{Cix) is continuous 
on the surface Ei = {x S | CiX — 0}. Indeed the jump of the vector field is equal to 2AiX = 
on Si. Moreover it is Lipschitz continuous as it is piececewise linear. Hence Proposition [1] applies. 

□ 

Corollary 4 Consider the differential inclusion in ^60]) . Suppose that (ii) and (iii) of Assump- 
tionUj) hold. Suppose that the multivalued mappings x g,i{x)sgn{hi{x)), 1 ^ i ^ m, are 
hypomonotone. Then for any initial data the differential inclusion 1160]} has a unique solution 
X : (0, T) IR^ that is Lipschitz continuous. 

The proof is straightforward and is omitted. 

Example 8 The mapping F : IR ^ M, x ^ sgn{aictaii{x)) , is hypomonotone with any 
k > -2- 

Remark 2 As noted in f55^ chattering may be due in sliding mode control appplications to the 
presence of parasitic dynamics. Simple modelling of these yield when inserted in (Qp the differential 
inclusion (see J551 (7) (8)]) 

\ / x{t) \ ( : 

1 I I a;,,(<) - I sgn{Cx{t)) (61) 

Xs{t) / V 
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with C = (0 1 0). The relative degree of the triplet {A,B,C) of this system is r = 3, where 

/ \ 

B = (1 0)^ and A = I 1 . This system does not fit within the above classes 
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of inclusions. Similar conclusions hold for the other form of parasitics in f55\, (3) (4)]- Such 
parasitics may be seen as a non collocation issue, that is known to greatly influence the stability 
of systems and usually may yield instability. The mere existence and uniqueness of solutions 
for such relative degree 3 systems is not trivial. In 142^ an example is given that possesses an 
infinity of absolutely continuous Filippov's solutions, but a unique so-called forward solution. One 
interesting question is to determine what kind of solution is approximated by the backward Euler 
method applied to f61]) which, according to 14 2[ Theorem 1] has a unique forward solution since its 
leading Markov parameter is CA^B ^ ^ > 0. A possible solution for this non collocation issue is 
the observer design of Example^ 

The differential inclusions in ([59|) - ((60|) are therefore discretized as follows: 

e f{tk,Xk) - pxk - Y.T=ii^i^k+i + B,)sgn{CiXk+i + D,) + p,Xk+i, a.e. on (0,T) 

(62) 

^''+''~^'' e f{tk,Xk) - pxk - g{xk+i)Sgn{h{xk+i) + pxk+i. a.e. on (0, T) 

(63) 

x(0) = Xo 

where the pi are the hypomonotonicity constants and X^Zli Pi — P- The result of Proposition [2] 
applies to ^ and feS]) . 




7.1 A simple hypomonotone case 

As shown in Section[2]on a simple monotone example, in practice the intersection may be computed 
as follows. Let us now illustrate this on the following system with hypomonotone multivalued part: 

x{t) e -{x{t) + l)sgn(.T(t)) + u{t) (64) 

with x{t) G IR. We may discretize it as: 

-{xk+i -Xk- huk - hpxk) e h{xk+i + l)sgn(a:fc+i) + hpxk+i, k^O, xq ^ x{0) (65) 

Notice that we may rewrite ((65|) as 

e Xk+i - {xk + huk + hpxk) + h{xk+i + l)sgn(xfc+i) + hpxk+i (66) 

Let us denote the mapping in the right-hand-side of l(66|) as F{xk+i). The set- valued mapping 
F{-) is strongly monotone \24\ Definition 2.3.1] for all p ^ 1. It follows from |24[ Theorem 2.3.3] 
that the generalized equation G F{xk+i) has a unique solution. 

For u{t) = 0, the following Lemma extends the Lemma [TJ 

Lemma 5 For alll > h > andxo G M, there exists fco such that Xko+n = and ''°+"+^^ feo+n 

for all n ^ 1. 

Proof:If xq G [—h,h], then fco = 0. Otherwise, for k < kg and Xk 4- ["^I'l]) the solution is 
given by : 

Xk - /isgn(xfc) , , 

1 + ft.sgn(a;fc) 
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From the solution ((67l) . the step fco for which Xk„ G [—h,h] can be easily found. Let us now 
consider that Xkg G [— /i, h]. The only possible solution for 

Xko + l - Xko = -H^ko + l + l)sfco + l 

€ sgn(a;fco+i) 

is Xko+i ~ and Sfc„+i ~ ~J^' next iteration, we have to solve 

Xko+2 = —hSko+2 

(69) 

and we obtain XkQ+2 ~ and Sko+2 = 0. The same holds for all Xko+n,Sko+n, n ^ 3, redoing the 

same reasoning. Clearly then the terms '^0+"+^ feo+n ^ppj-Qxiniating the derivative, are zero 

h 

for any h > 0. □ 
We conclude that in this case also the system and its derivative are correctly approximated at 
the zero value on the sliding surface. There is no spurious oscillation around the switching surface. 



8 Detailed Implementation of Implicit Euler Discretization 
of the general case (1601) 

This section is devoted to the implementation and the study of the numerical algorithm. The 
interval of integration is [0,T], T > 0, and a grid to = 0, t^+i = tk + h, k ^ 0, = T is 
constructed. The approximation of a function /(•) on [0,T] is denoted as f^{-), and is a piecewise 
constant function, constant on the intervals [tk,tk+i). We denote f^{tk) as fk- The time-step is 
/i > 0. 



8.1 Time— discretization 



x{t)^ f{x{t),t)-g{x{t))s{t) 
y{t) = h{x{t)) 
s{t) e Sgn(y(0) 



(70) 



Starting from l|60p . let us introduce a new notation. 



where s[t) £ and y{t) £ are complementary variables related through the Sgn(-) multi- 
valued mapping. According to the class of systems lfT3|). l(59|) or (|60|) that we are studying the 
functions /(•) and g{-) are defined either in a fully nonlinear framework or by affine functions. We 
present the time-discretization in its full generality and specialize the algorithms for each case in 
Section [Ql 

Let us now proceed with the time discretization of ((70l) by a fully implicit scheme : 

Xk+i = Xk + hf{xk+e,tk+e) - hg{xk+j)sk+i 

yk+i^h{xk+i) (71) 
Sk+i £ Sgn{yk+i) 

where Xk+e = 9xk+i + (1 - 0)xk, Xk+^ = ^Xk+i + (1 - l)xk, and tk+e = Otk+i + (1 - 0)tk+i, with 
9 = [0, 1] and 7 £ [0, 1]. As in [l], we call the problem ((7T|) the "one-step nonsmooth problem". 

This time-discretization is slightly more general than a standard implicit Euler scheme. The 
main discrepancy lies in the choice of a ^-method to integrate the nonlinear term. For = 0, we 
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retrieve the explicit integration of the smooth and single valued term /(•). Moreover for 7 = 0, the 
term <?(•) is explicitly evaluated. The flexibility in the choice of 9 and 7 allows the user to improve 
and control the accuracy, the stability and the numerical damping of the proposed method. For 
instance, if the smooth dynamics given by /(•) is stiff, or if we have to use large step sizes for 
practical reasons, the choice of 6* > 1/2 offers better stability properties with respect to h. 



8.2 Mixed Complementarity Problem 

The so-called "one-step nonsmooth problem" appears to be a Mixed Complementarity Prob- 
lem (MCP) that we have to solve at each time-step. Let us deflne what is a MCP : 

Definition 2 (MCP) Given a function f : M"^ IR'' and lower and upper hounds l,u d M'^ , 
find z e TR? , w, tj G IB?^ such that 

F{z) = w — V 

I Z ^ U 

{z~l)^w^O 
[u — z)'^v ~ 

where M = MU {+00, —00}. 

Note that the problem l(72|) implies that 

-F(z)e7V[,,„](z). 



(72) 



(73) 



The relation l(73|) is equivalent to the MCP l(72|) if we assume that w is the positive part of F{z), 
that is w = F^{z) = max{0, F{z)) and v is the negative part of F{z), that is t; = F~{z) — 
max(0, —F{z)). 



The One step nonsmooth problem as a MCP Let us deflne the MCP by 



F{z) = 



Xk+l 
Sfc+l 



Xk+l - hf{xk+e,tk+e) + hg{xk+^)sk+i 
~h{xk+i) 



(74) 



—00, i = 1 . . .n 
— l,i = n+ l...TO 



-00, i = 1 . . .n 
- 1 , z = n + 1 . . . TO 



If z solves the MCP (|74|) . the bounds u and I and the condition (z — lY'w = 0, (u— z)'^v — imply 
that 



0, i = 1 ... 71 
y~ 0,i — n + 1 . . .m 

The MCP is then given by 

Xk+l - Xk- hf[xk+e) + hg{xk+^)sk+i 
~h{xk+i) 



Q,i = l...n 

^ 0,i = n + 1 . . .m 





y~{xk+i) - y+{xk+i) = -y{xk+i) 



(75) 



-1 < Sk+l < 1 

(sfe+i 4- l)'^y~{xk+i) = 
{l-.sk+ify+{xk+i) = 



(76) 



It is clear that the problem ((7T|) is equivalent to the MCP defined by l(7i|) . The results of existence 
and uniqueness of solution of ([72l) or equivalently l(73|) are related to the monotonicity properties 
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of F (P-properties or coherent orientations of the associated afHne map for particular cases of 
bounds and afHne function F(-). Without entering into further details, we refer to [30^ [24] for the 
main results. 



Numerical Solvers The MCP l(72|) can be solved by a large family of solvers based on Newton- 
type Methods and interior-points techniques. We refer to [TT] for a comparison of several solvers 
based on Newton's method. The numerical implementation of the MCP solvers are often based 
on the computation of the Jacobian matrix of the function F(-) with respect to z. The Jacobian 
matrix is explicitly given in the case defined in (fM)) by 



V,F(z) 



/ - h0Vxfixk+e, tk+e) + hjVxg{xk+'r) ® Sfc+i hg{xk+^) 
Vxh{xk+i) 



(77) 



where ® denotes the simple contracted tensor product and the third-order tensor xg{x) is the 
Jacobian of g with the respect x given by the following component: 

/Y7 I \\ dgki{x) 

{^xg{x))kip = —g^ — ■ (78) 

8.3 Newton's linearization and Mixed Linear Complementarity Prob- 
lems 

Due to the fact that two of the systems that are studied in this paper involve affine functions /(•) 
and g{-), we propose to solve the "one-step nonsmooth problem" (fzT]) by performing an external 
Newton linearization, which yields a Mixed Linear Complementarity Problems (MLCP). 

Newton's linearization The first line of the problem l|7ip can be written under the form of a 
residue TZ depending only on Xk+i and Sk+i such that 

7^(xfe+l,sfc+l) = (79) 

with 7l{x, s) = x — Xk — hf{0x+ (1 — 0)xk, tk+e) + hg{'-fx+ (1 — j)xk)s. The solution of this system 
of nonhnear equations is sought as a Hmit of the sequence {x"^^-^, Sfe+i}Q,gJV ^^^^ ^^^^ 

Xk+i ~ Xk 

nUxttl^tl) = ^k+i - ^fc - hfix-,^,) + V,7^(x^+l, s^+i)(x^+i - x^+i) + hgixt+^)s'^,tl = « 

(80) 

The computation of the Jacobian of TZ with respect to x, denoted by M{x, s) leads to 

M{x, s) = V^7^(x, s) = I- heVxf{Ox + (1 - 6)xk,tk+e) + h-fV xg{-ix + (1 - 7)xfe) ® s. (81) 
At each time-step, we have to solve the following linearized problem, 

x^+i - Xk - hf{xt_,e,tk+e) + A/(x^+i, s^+l)(^fe^^ - ^fe+i) + hg{xt+^)slXl = 0, (82) 

that is 

^IXI = + M~\xt+,,st+,) [xk - xt+, + hfixt+g,tk+e) ~ hg{xt+^)s^ktl] ■ (83) 

The matrix M is clearly non singular for small h. The same operation is performed with the 
second equation of l|7ip leading to the following linearized equation 

y^tl = Vk+i + V./i(x^+i) [xlXl ~ xl_,^] (84) 
Inserting ijSS]). we get the following Hnear relation between y'^Xi ^'^^ */e+i ' 

vtXl = Vk+i + ^M^Ui) \M-\x'^,^„st+^)ixk - x^+i + hfix'^,^„tk+e) - hg{xl^^)slXl) 
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Mixed linear complementarity problem (MLCP) To summarize, the problem to be solved 
in each Newton iteration is: 



s 



(86) 



with G K"'''" and b e JR" defined by 

(87) 

KXl = vUi + ^Mxt+i) [M-\x'^,+„.st+,)ixk - x^+i + hf{xt+e,tk+e))] 

The problem ((86l) is equivalent to a MLCP which can be solved under suitable assumptions 
by many linear complementarity solvers such as pivoting techniques, interior point techniques and 
splitting/projection strategies. The reformulation into a standard MLCP follows the same Hne as 
for the MCP in the previous section. One obtains, 

yttl'^ vttl = -w:tl4tl + Ktl 

As for the MCP, there exists numerous methods to numerically solve MLCP. In the worst case 
when the matrix W^X^X ^as no special properties, the MCLP can be always solved by enumerative 
solvers for which various implementations can be found. With some positiveness properties [2i] . 
standard methods for LCP[19] can be straightforwardly extended. Among these methods, we can 
cite the family of projection/splitting methods, interior point methods and semi-smooth Newton 
methods (see [l] for an overview). 

8.4 The special cases of the afRne systems 

In this section, we specify the time-discretization to the two other classes of systems ifTS]) and 
((59|) and for particular value of 9 and 7. 



8.4.1 Time— discretization of the system (I13l> 

For the system (fT3)) . the function g{x) is reduced to the matrix B and the function h{x) is afHne, 
that is h{x) = Cx + D. The matrix W^X^X and b'XX\ are then given by 

WXtX = hCM-\xt+„st+„tk+i)B 

(89) 

Ktl = Vt+i + ^ [M-'{xt+iM+i){xk - xl^, + hf{xl^,,tk+e))\ 
with 

M{x,t) = I - hev ^f{ex + {I ^ e)xk,et + {\ - e)tk) (90) 

\{ CB > Q then the matrix W^X^X ^^^^ positive definite for sufficiently small h. This result 
ensures the existence and uniqueness of a solution of the MCLP. Furthermore, standard pivoting 
techniques such as Lemke's method or projection/splitting such as Projected Successive Over- 
Relation (PSOR) compute the solution. 

Semi-implicit discretization with = The matrix M is then reduced to 

M{x,t)=I (91) 
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and 

W = hCB 

(92) 

hk+i = yk + hCf{xk,tk) 

In this particular case, there is no need to perform some Newton iterations because the system to 
be solved at each time-step is linear. Furthermore, the MLCP has a solution for any h > under 
the assumptions that CB > 0. 

Fully implicit discretization with an afflne function The same conclusion can be drawn 
if VxfiOx + (1 — 9)xk, 9t + (1 — 6)tk) is equal to a constant matrix E that is when /(•) is linear 
time-invariant and given by f{x,t) — Ex{t) + F. In this case, the matrix M reduces to 

M{x, t)=I- hOE (93) 

and 



W = hC{I - heE)-^B 
bk+i = Vk + hC[I - hOE)-^ [Exk + a] 

8.4.2 Time discretization of the system (15911 

For the system l[59|) . we recall that g{ ) and h{-) are given by 

g{x) = [g,j {x) = Ajx + Bj,j ^l...m]e JR" 
h{x) = [hi{x) = CiX + D,,i^l...m](^IBJ^. 

The components of g(-) can be explicitly expressed by 



(94) 



(95) 



gu {x) ^^AukpXp + Bi, (96) 

and therefore, the Jacobian of g{-) is given by 

, , dgkijx) 

{Vxg{x))kip = —g^ — = Ai^kp- (97) 

The Jacobian of h takes the following simple form: 

Vh{x) = C=[C,,2 = l...m]eJR'"^". (98) 
After the newton linearization, we have to solve at each Newton's loop the MLCP l|86p with 

W^+,' = hCM-\xt+,,st+„tk+i)gixt+^) 

(99) 

Ktl = Vt+i + C [M-\xl^,,tk+i){xk - x^+i + hf{xl^,,tk+e))\ 
Semi— implicit discretization with = 7 = If 7 and 9 vanish, the residue Ti given by 

Tl{x, s) = X - Xk - hf{xk,tk) + hg{xk)s (100) 

is linear in x and s. In this particular case, there is no need to perform Newton's iterations. The 
MCLP defined by l(99| can be simplified to 

Wk+i=hCg{xk)bk+i=Vk + C[hf{xk,tk))] (101) 
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Algorithm 1 Implicit Euler time-discretization with a generic MCP solver 
Require: System definition: f, g, h 
Require: x(0) the initial condition 
Require: to, T time-integration interval 
Require: h time-step 

Require: 0,7 numerical integration parameters 
Ensure: ({xj, {sj, {yk}), k e {1, 2, . . .} 



k^O; xo^x(O); yo ^ y(0) = h(x(0)); tauo ^ 
// Time integration loop 
while tk < T do 

Solve the MCP ([76]) for x^+i, s^+i, yk+i with F, I and u given by ([74]) and the Jacobian \/zF{z) 

given by ([TTll. 

//Update 

Xk ^ Xk+i; Sk ^ Sk+i; yk ^ yk+i 
/ / time iteration 
tk^tk+i; k^k+1 
end while 



Algorithm 2 Implicit Euler time-discretization with an external Newton loop and a MLCP solver 
Require: System definition: f, g, h 
Require: x(0) the initial condition 
Require: to, T time-integration interval 
Require: h time-step 

Require: 6,^ numerical integration parameters 
Require: e Newton's method tolerance 
Ensure: ({xk}, {sk}, {yk}), k e {1, 2, . . .} 



k^O; xo^x(O); yo ^ y(0) = h(x(0)); tauQ ^ 

// Time integration loop 
while tk < T do 

(T ^ 0- x° ^ x°- s° ^ s°- v° ^ v° 

/ / Newton s loop 

while ||7^(x^+l,s^+l)|| >£do 

M-Hx^+i,s^+i) ^ (I - h0V,f(x«^_„tk+i) - h7V,g(x«+,)®s^+i)-i. 

Kti ^ Vk+i + Vxh(x^+i)M-Hs^+i) [xk - x^+i + hf(x^_,J] 
Solve the MLCP ^ for yki\sk+i^ 

^ \+l ' Yk+l ^ X|^_|_^ 

^ x^+1 + M-i(s^+i) [xk - x^+i + hf(x?+i) -I- hg(x^+i)s^+i] 
a ^ a + 1 
end while 
//Update 

„ ^ „a . e ^ „Q!+1. ,, . .,"+1 

Xk+1 ^ X|^_|_^, Sk+1 ^ s^^;^ , yk+1 ^ yk+1 
/ / time iteration 
tk^tk+i; k^k+1 
end while 
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8.5 Algorithms 

We propose in this section two algorithms to sum-up the numerical implementation of the implicit 
Euler time-stepping scheme. The Algorithm [T] describes the implementation with a generic MCP 
solver and the Algorithm [2] describes the numerical implementation of the algorithm with an 
external Newton linearization and a MCLP solver. 

In the case of the system (flSl) with an afHne function /(•) or 6* = 0, the MLCP matrix W can 
be computed before the beginning of the time loop, saving a lot of computing effort. In the case 
of the system ((59|) with 6* = 7 = 0, the MLCP matrix W can be computed before the beginning 
of the Newton loop. 
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9 Numerical experiments 

Let us illustrate the above developments with numerical integrations performed with the SICONOS 
software platform of the INRI^l [H [3] which is designed for the simulation of multivalued nons- 
mooth systems. 

9.1 Chattering free stabilization 

Let us consider the following continuous-time closed loop system from |26] given by 



" 1 ■ 




■ " 


-ci 


X — 


a 



sgn([ ci 1 ] x). 



(102) 



As it is shown in [26] the trajectories obtained by an explicit Euler discretization exhibit spurious 
oscillations which are described by period-2 cycle around the sHding manifold. On Figure [H the 
trajectories obtained by implicit discretization are shown using the implicit Euler discretization 
with h = 1, h = 0.3, h = 0.1 and h = 0.01 and with ci = 1 and a = 1. As it has been predicted 
by theoretical discussions of Section IH the sliding manifold is reached in finite time and without 
any chattering. Indeed, the matrix CB = a = 1 satisfies the assumptions of Lemma [3l Note that 
the algorithm is also very robust in the sense that the simulation can be performed with relatively 
large time-steps. 



9.2 Example [3]: Multiple sliding surfaces 

Let us consider the example [3l The system can be defined in the form ifTSj) with 



B = 



1 2 

2 -1 



C = 



1 2 

2 -1 



D = 0, f{x{t),t) = 



This example illustrates Lemma [3] since CB 



5 
5 



(103) 



The results displayed on Figure [5] show 



that the system reaches firstly the sHding surface 2x2 +xi ~ without any chattering, The system 
then sHdes on the surface up to reaching the second sHding surface 2xi — X2 — and comes to rest 
at the origin. 



9.3 Extensions to ZOH discretized systems 

The extension to ZOH discretized systems is illustrated on a first example taken from [57]. In the 
notation of Section [5^ the LTI system with an ECB-SMC controller is defined by the following 
data, 



F 





-0.1 



1 

-a2 



G 



C [ ci 1 ] 



(104) 



Starting from the initial data, xq — [0.55, 0, 55]-^, Galias and Yu [27] have shown that the Explicit 
ZOH discretization of the system with ai = —2, a2 = 2, ci = 1 and h = 0.3 exhibits a period-2 
orbit. The results are reproduced on Figure 6(a) On Figure [6 (b)[ the Implicit ZOH discretization 
as proposed in Section [5^ is free of chattering. On Figure [71 a comparison is given between the 
ZOH and the Euler discretization of the vector field /(•, •). Another example taken from [57] in 
the MIMO case is given by the following parameters. 











1 " 




" 


" 


F = 


1 


1 


1 


, G = 


1 







-1 


-3 


1 







1 



c 



1 1 
1 1 



(105) 



Similar results are depicted on Figure [8l 

^http : //siconos .gf orge . inria. f rZ] 
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xi(t) 



(a) h = 0.3. Explicit Euler 



(b) ft = 0.1. Explicit Euler 
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xi(t) 



0.5 0.6 0.7 



(c) h = 1. Implicit Euler 
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"id) 

(d) h = 0.3. Implicit Euler 
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(e) h = 0.1. Implicit Euler 
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(f) h = 0.05. Implicit Euler 



Figure 4: Equivalent control based SMC, ci = 1,Q! = 1 and xq = [0,2.21]^. State xi{t) versus 

X2it). 



















































\ 

\ 



(a) state xi{t) and X2(t) versus time (b) phase portrait X2{t) versus xi{t) (c) sgn function si{t) and S2{t) 
Figure 5: Multiple Sliding surface, h = 0.02, a;(0) = [1.0, -1.0]'^ 
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-0.1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.1 0.2 0.3 0.4 0.5 0.6 0.7 



xi(t) xi(t) 

(a) h = 0.3. Explicit ZOH (b) h = 0.3. Implicit ZOH 



Figure 6: Equivalent control based SMC, ai ~ —2, a2 = 2, ci = 1 and h = 0.3. .tq = [0.55, 0, 55]"^ 
State xi[t) versus X2{t). 
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(a) Explicit implementation 
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-0.1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

xi(t) 

(b) Implicit implementation 



Figure 7: Comparison of Euler and ZOH discretizations of ECB-SMC system, ai — —2, 02 — 2, 
ci = 1 and h = 0.3. = [0.55,0,55]"^ State xi{t) versus X2{t). 




(a) h = 0.3. Explicit ZOH (b) h = 0.3. Implicit ZOH 



Figure 8: Equivalent control based SMC, Trajectory with initial point xq = [0.05,-0.5,0.02 
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9.4 Lyapunov-based robust control 

We propose in this section to give an numerical example which fits with the example (fT8| of a 
Lyapunov-based discontinuous robust control. Let us consider the following system 



b{t) = -x{t) ^ u{t) + j{t) 



(106) 



with j{t) ~ asin(t) and u{t) — sgn(a;(t)). Is is obvious that, as expected, the impHcit method 
yields a smooth stabilization at a; = whereas the expHcit Euler has significant chattering. Fig- 
ure 9(c) illustrates the fact that the controller varies inside the multivalued part of the sgn function 



in order to assure the existence of an equilibrium point. 



(a) State xi{t) vs. time, h = 0.1. Implicit (b) State xi(t) vs. time, h = 0.1. Explicit 
Euler Euler 




time t time t 

(c) Control u{t) vs. time, h = 0.1. Implicit (d) Control u{t) vs. time, h = 0.1. Explicit 
Euler Euler 



Figure 9: Lyapunov-based discontinuous robust control, h = O.l a = 0.1 



9.5 The Filippov example 

Example 9 Let us consider now the well known Filippov example which can be defined in the 
form (03) with 

D = 0, f{xit),t)=0. (107) 

The trajectories may slide on the codimension 2 surface given by Cx = 0, that x ~ 0. 

The results displayed on Figure 19.51 show that the system reaches the origin without any 
chattering. The sufficient conditions of the Lemma [3] are not satisfied when seems to indicated 
that these conditions has to be improved. 

9.6 Example [H Observer based SMC 

Let us illustrate the performance of our implementation on the observer based SMC described by 
the Example m The dynamics is given by lfT7|) with fc = 1 and r = 0.001. The initial conditions 



B 



1 
1 
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(a) state x-i(t) and X2{t) versus time 
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0.2 0.4 0.6 0.8 



(b) phase portrait X2{t) versus x\{t) 



(c) sgn function s\{t) and S2(t) 

Figure 10: Multiple Sliding surface. Filippov Example, h = 0.002, x{Q) = [1.0,-1.0]'^ 
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are chosen as [2.0,0,0,0]^. The numerical parameters are given by h ~ 0.1 that is a sampling of 
lOHz and 9 = ^ 1. On Figure [TT| the error between the reference command and the observer 
state is given. On Figure [T2l we can observe the behavior of the command without any chattering. 

0.0035 
0.003 
0.0025 
0.002 
0.0015 
i 0.001 
0.0005 


-0.0005 
-0.001 
-0.0015 

2 4 6 8 10 

time 

Figure 11: Observer based SMC: Error e{t). fc = 1 and t = 0.001. /i = 0.1 6* = 1, 7 = 1 
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I 0.4 

o 

0.2 





-0.2 

02468 10 02468 10 

(a) control u{t) = Sgn{Cx{t)) (b) Zoom on control u{t) = Sgn{Cx(t)) 

Figure 12: Observer based SMC: Control. A: = 1 and r = 0.001. h = 0.1 6 = ^ 1 

On the three previous examples, one sees a very accurate and smooth stabilization on the 
sliding surface, even for values of h not so small. 

Influence of the integration parameters 6 and 7 In the following numerical experiments, 
we discuss the role of the numerical parameters and 7. Due to the fact that the function 
g{-) is linear and reduced to a matrix B = [1,0,0,0]^, the parameter 7 has no influence on 
the numerical time-integration. On the contrary, the parameter 9 has a huge influence on the 
stability of the integration. Indeed, the implicit Euler integration (9 = 1) of the smooth term 
/(•) is unconditionally stable. This is not the case for the explicit Euler 9 — and for the 
chosen parameters k and t, the instability of the scheme for ft, = 0.1 does not allow to proceed to 
integration. On Figure [131 the instability of the scheme is illustrated and appears as a chattering 
on the state x. The stability is retrieved for h < 0.005. 

Other choices of 9 can be made to improve the numerical time-integration of the smooth dy- 
namics. For instance, 9 = 1/2 yields a second order scheme for the integration of /. Unfortunately, 
the scheme is not of second order due to the fully implicit integration of the multi- valued part. 
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Nevertheless, it can be interested to use of 6* = 1/2 to decrease the numerical damping of the 
scheme on the smooth term. On the Figure [TH two simulations of the observer based SMC are 
presented for = 1 and = 1/2. The parameter r has been modified to 0.5 to correctly integrate 
the parasitic dynamics with the same time-step h = 0.1. 




(a) state xi{t) (b) control u{t) = Sgn{Cx{t)) 

Figure 14: Observer based SMC with 6* = 1 and 61 = 1/2. fc 1 and r 0.5. h = 0.01,7 1 

10 Conclusions 

In this paper the backward Euler method is studied on specific classes of Filippov's systems 
that encompass sliding mode control systems. It is shown that such implict schemes allow a 
smooth accurate stabilization on the sliding surface, even in case of codimension larger than one. 
Despite the backward Euler method has been studied and used for a long time in other fields 
like contact mechanics and electric circuits simulation [T], it seems it has not yet been used in 
the sliding mode control community. This work therefore constitutes the introduction of a new 
discretization method for EBC-SMC systems. The novelty compared to numerical simulation is 
that this time one has to consider not only the numerical simulation, but also the implementation 
on real processes. Perhaps one obstacle to the dissemination of the method is that at first sight, 
the controller designed from a backward philosophy looks like a non causal controller. However 
as shown in this paper this is not the case. This paper paves the way towards the study of a new 
family of discrete-time sliding mode controllers. 
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